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TIME  DEPENDENT,  COMPRESSIBLE  SIMULATIONS  OF  SHEAR  FLOWS: 

TESTS  OF  OUTFLOW  BOUNDARY  CONDITIONS 

I.  INTRODUCTION 

this  paper  describes  recent  numerical  simulations  performed  at  NRL  of 
shear  flows  transitioning  to  turbulence,  and,  in  particular,  the  development 
and  evolution  of  coherent  structures.  There  are  two  aspects  of  this  problem 
which  are  addressed  in  this  report.  The  first  is  developing  and  testing  the 
model  that  was  used  in  these  studies.  In  particular,  we  are  concerned  with 
the  treatment  of  inflow  and  outflow  boundary  conditions  suitable  for  both 
compressible  and  incompressible  flows .  The  second  aspect  is  using  thi3  model 
to  describe  shear  flows. 

The  first  problem,  developing  the  proper  computational  tools,  has  been 
the  major  goal  of  the  last  year's  work  and  is  discussed  in  Sections  III 
through  VI.  The  numerical  model  we  are  usinq  now  is  a  restructured  version 
of  the  FAST2D  computer  code.  This  incorporates  the  Flux-Corrected  Transport 
( FCT)  continuity  equation  algorithm  (Boris,  1976b;  Boris  and  Book,  1976) 
which  has  been  tested  extensively  for  shock,  detonation,  and  beam-qenerated 
turbuelnce  calculations  (e.g.,  Book  et  al.,  1980,  Oran  et  al.,  1982;  Picone 
and  Boris,  1983).  Since  the  algorithm  is  explicit  in  its  present  form,  the 
code  is  particularly  efficient  for  studying  flows  that  move  at  a  substantial 
fraction  of  the  speed  of  sound  in  the  material.  Using  time  step  splitting 
techniques,  we  couple  FCT  to  alcorithms  for  the  other  physical  processes  we 
want  to  represent.  In  Section  III  below  we  describe  the  general  features  of 
the  code.  In  Section  IV,  V,  and  VI  we  describe  tests  of  the  new  outflow 
boundary  conditions. 

The  second  aspect,  application  of  the  model,  is  the  goal  of  the  upcominq 

research.  We  have  used  the  model  to  simulate  time  dependent  flows  in  two 

confiqurations,  the  splitter  plate  and  the  round  jet,  for  which  substantial 
Manuscript  approved  October  28,  1983. 
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data  exist  on  the  transition  to  turbulence.  At  the  end  of  this  paper,  in 
Section  VII,  we  describe  some  preliminary  results  and  their  implications. 
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II.  BACKGROUND 


The  development  and  structure  of  turbulent  flows  is  the  focus  of  intense 
study,  we  new  Know  that  flows  which  previously  were  thought  to  he  totally 
chaotic  and  statistical  in  nature  are  dominated  by  the  persistence  of  rela¬ 
tively  large  structures.  These  coherent  structures  were  previously  ignored 
since  their  existence  was  masked  or  de-emphasized  by  experimental  averaging 
techniques.  The  classical  description  of  turbulence  and  the  mechanisms 
responsible  for  its  development  are  now  considered  deficient  in  their  expla¬ 
nations  of  these  transient  but  organized  and  persistent  structures. 

The  classical  description  of  turbulence  evolved  from  the  observed 
behavior  of  fluid  flows  as  a  function  of  Reynolds  number.  Many  flows  exhibit 
a  series  of  sudden  transitions  to  new  flow  patterns  as  the  Reynolds  number  is 
increased.  Each  transition  results  in  an  increasingly  complicated  flow,  and 
at  high  Reynolds  number  the  flows  become  irregular  and  appear  chaotic  in  both 
space  and  time.  The  transitions  to  the  succession  of  flow  patterns  may  be 
caused  by  a  seouence  of  fluid  instabi lities ,  each  of  which  breaks  some  sym¬ 
metry  in  the  previously  stable  flew  pattern  and  introduces  some  new  scale  in 
the  flow  pattern  (Liepmann,  1979) .  The  transition  to  a  purely  chaotic  flow 
was  postulated  to  occur  through  an  infinite  succession  of  instabilities,  each 
contributing  to  the  increasing  frequency  content  of  the  flow  (Landau  and 
Lifshitz,  1959). 

Several  experimental  observations  have  seriously  eroded  confidence  in 
the  completeness  of  such  a  turbulence  model  for  physical  flows.  First, 
although  laboratory  created  grid  turbulence  comes  close,  no  flow  has  yet  been 
shown  to  exhibit  pure,  homogeneous  isotropic  turbulence  in  the  classical 
sense.  Such  a  state  is  really  a  limiting  condition.  Further,  since  coherent 
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structures  exist  on  the  larger  scales,  it  is  reasonable  to  assume  the  exis¬ 
tence  of  similar  structure  on  all  scales  laraer  than  the  dissipation  scale. 
Second,  the  transition  to  turbulence  does  not  occur  through  am  infinite  suc¬ 
cession  of  instabilities,  but  after  the  appearance  of  relatively  few  insta¬ 
bilities,  typically  three  or  four.  This  has  led  to  the  postulate  of  a  theo¬ 
retical  connection  between  strange  attractor  theory  and  transition  to  turbu¬ 
lence.  Third,  intermittency  in  free  turbulent  shear  layers  indicates  the 
presence  of  a  thin,  sharp  interface  between  turbulent  fluid  and  irrotational 
fluid.  This  finding  leads  to  questions  about  whether  such  sharp  interfaces 
can  be  represented  as  a  diffusive  effect.  Finally,  the  recent  discovery  that 
coherent  structures  dominate  flows  which  were  previously  believed  to  approach 
pure  turbulence  directly  causes  us  to  reconsider  the  basic  assumptions  in  the 
classical  theory  of  turbulence. 

The  search  for  new  concepts  of  the  nature  of  turbulence  has  centered  on 
understanding  the  mechanisms  of  the  transition  to  turbulence  in  several 
simple  fluid  flows.  Shear  flows  generated  by  splitter  plates  and  round  jets 
(e.g..  Brown  and  Roshko,  1974;  Browand  and  Weidman,  1976;  Roshko,  1976; 
Fiedler  and  Wygnanski,  1970)  exhibit  all  of  the  troublesome  intricacies  asso¬ 
ciated  with  the  transition  and  in  addition  are  particularly  appropriate  for 
practical  applications.  Both  larger  and  smaller  scales  of  chaotic  fluid 
motion  develop  in  these  systems.  Flow  near  the  orifice  may  initially  be 
laminar  for  even  high  Reynolds  numbers,  and  several  distinct  transitions  are 
easily  discerned  before  the  onset  of  an  apparently  fully  turbulent  shear 
layer.  It  is  still  quite  difficult,  however,  to  pinpoint  at  which  point  in 
the  flow  the  label  "turbulent"  is  applicable.  Indeed,  coherent  structures  of 
large  scale  have  been  found  to  dominate  well  downstream  in  flows  which  appear 
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chaotic  on  smaller  scales.  By  using  similarity  arguments,  for  example,  we 
see  that  the  splitter  plate  flow  is  always  dominated  by  ever  larger  coherent 
structures.  Through  flow  visualization  experiments,  a  great  deal  of  insight 
has  been  achieved  into  the  exact  mechanisms  involved  in  the  individual  insta¬ 
bilities  as  well  as  some  indication  of  the  sequence  of  appearance  of  various 
scales  of  motion. 

Although  the  flows  near  the  orifice  may  be  laminar  for  both  jets  and 
splitter  plates,  the  shear  layer  generated  is  unstable  to  the  Kelvin- 
Helmholtz  instability.  Small  perturbations  in  the  flow  grow  into  nonlinear 
waves  which  break  and  roll  up,  transforming  the  original  vorticity  of  the 
shear  layer  into  isolated  clumps.  The  primary  wavelength  generated  by  the 
instability  is  usually  that  of  the  fastest  growing  mode  for  that  particular 
geometry  or  of  some  impressed  wavelength  determined  by  boundary  conditions 
or  initial  conditions.  Further  development  of  the  shear  layer  proceeds 
through  the  pairing  of  vortex  clumps,  a  process  which  may  be  repeated  many 
times  downstream.  One  of  the  effects  of  pairing  is  to  generate  subharmonics 
of  the  original  unstable  wavelength,  but  smaller  wavelength  disturbances  are 
created  as  well.  These  disparate  wavelengths  arise  from  at  least  two  causes: 
imperfect  pairing  due  to  small  fluctuations  in  the  flow  leaving  an  unpaired 
vortex  which  then  merges  with  a  previously  formed  pair  (Browand  and  Winant, 
1973),  and  the  generation  of  small-scale  disturbances  in  the  interaction  of 
the  cores  of  the  two  vorticity  clumps  (Zabusky,  1981;  Overman  and  Zabusky, 

1981 ) . 

It  should  be  emphasized  that  the  flow  may  remain  two-dimensional 
throughout  this  process,  and  that  the  rich  frequency  content  of  the  flow  has 
been  produced  by  the  action  of  only  a  few  two-dimensional  instabilities.  As 


a  matter  of  practical  importance,  however,  the  flow  almost  always  becomes 
three-dimensional  when  the  Reynolds  number  is  high  enough.  Three  dimensional 
instabilities  generally  begin  to  play  a  role  about  the  same  time  that  visible 
vortex  pairing  ceases.  This  is  primarily  because  the  two-dimensional  vortex 
rolls  or  loops  are  deformed  by  wave  instabilities  in  the  thiru  dimension  or 
by  the  presence  of  boundary  conditions  imposed  by  the  system  size.  As  these 
deformations  grow,  further  two-dimensional  pairing  becomes  difficult  to 
observe  and  growth  in  the  third  dimension  seems  to  be  preferred.  The  non¬ 
linear  growth  phase  for  the  three-dimensional  instabilities  is  marked  by 
increasing  entanglement  of  vortex  lines  or  loops,  and  the  flow  becomes  more 
strongly  turbulent. 

The  shear  layer  grows  as  vortex  cores  entrain  irrotational  fluid  into 
the  shear  region.  On  the  macroscopic  scale  the  random  walk  of  long  vortex 
cores  or  filaments  spreads  the  reqion  containing  the  vorticity.  This  process 
resembles  an  eddy  diffusion  of  vorticity  at  scales  longer  than  the  dominant 
scale  of  the  eddies  responsible  for  the  transport.  Representing  this  convec¬ 
tive  phenomenon  as  microscopic  diffusion  at  short  scales  is  incorrect,  how¬ 
ever,  because  random  walk  mixing  increases  these  gradients  rather  than 
reducing  them.  Hie  eddy  diffusion-mixing  length  models  work  surprisingly 
well  macroscopically  because  the  diffusion  is  nonlinear,  and  the  eddy  coef¬ 
ficient  is  large  where  the  macroscopically  averaged  vorticity  is  large.  The 
front  which  propagates  from  such  nonlinear  diffusion  can  be  sharp  and  hence 
has  the  potential  to  represent  the  thin  transition  between  rotational  and  ir¬ 
rotational  fluid.  True  vorticity  diffusion  enters  into  the  picture  only  at 
the  smaller  scales,  transporting  vorticity  into  the  continually  thinning 
irrotational  layers  of  fluid  which  have  been  engulfed. 
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In  the  work  presented  below,  we  describe  numerical  simulations  of  two- 
dimensional  shear  flows  which  go  unstable  through  the  Kelvin-Helmholtz 
instability.  Our  coal  is  to  produce  accurate  enough  calculations  of  the 
time-evolution  of  the  major  physical  Quantities  so  that  we  can  analyze  the 
flow  behavior  and  test  our  concepts  of  the  transition  to  turbulence.  We 
consider  two  geometries:  cylindrical,  such  as  the  round  jet,  and  Cartesian, 
such  as  the  splitter  plate.  In  particular,  we  are  interested  in  the  initial 
transient  behaviour  and  in  the  resulting  unsteady  pattern  of  structures. 
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III.  THE  NUMERICAL  MODEL 

The  code  FAST2D  was  used  to  perforin  the  shear  flow  calculations 
described  below.  This  code  consists  of  a  solutions  of  the  time -dependent 
conservation  equations  for  mass,  momentum  and  energy  coupled  to  algorithms 
describing  gravity,  molecular  and  thermal  diffusion,  and  chemical  reactions 
with  energy  release.  Ttiese  various  parts  of  the  code  are  coupled  by  timestep 
splitting  and  can  be  turned  on  or  off  independently  by  logical  controls,  as 
required  by  the  problem  to  be  studied. 

The  continuity  equations  are  solved  using  the  FCT  algorithm  JPBFCT,  an 
advanced  version  of  ETBFCT  (Boris,  1976a).  FCT  is  a  finite-difference  tech¬ 
nique  for  solving  the  convective  equations  which  is  particularly  useful  in 
problems  where  sharp  discontinuities  arise  and  are  maintained  throughout  the 
calculations.  These  discontinuities  may  be  shocks  or  contact  surfaces.  In 
the  case  studied  here,  we  are  concerned  with  interfaces  in  material  density. 
The  algorithm  modifies  the  linear  properties  of  a  high  order  algorithm  by 
adding  sufficient  diffusion  during  convective  transport  to  prevent  dispersive 
ripples  from  arising,  and  ensures  that  all  conserved  quantities  remain  mono¬ 
tonic  and  positive.  This  added  diffusion  is  subtracted  out  appropriately 
where  not  needed  in  an  antidiffusion  phase  of  the  timestep  to  maintain  second 
order  accuracy. 

Another  important  feature  of  the  FCT  algorithms  is  their  ability  to 
divorce  the  arid  motion  from  the  fluid  flow.  This  freedom  has  been  used  to 
incorporate  variably  spaced  grids  as  well  as  adaptive  grids  (Book  et  al., 
1980;  Oran  et  al.,  1982)  which  automatically  follows  regions  where  there  are 
sharp  gradients  and  more  resolution  is  needed.  In  transition  to  turbulence 
calculations,  this  means  that  realistic  spatial  inflow  and  outflow  problems 
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can  be  solved  rather  than  the  more  idealized  periodic  problems  where  temporal 
but  not  spatial  problems  can  be  considered. 

The  JPBFCT  routine  solves  the  one-dimensional  continuity  equation  in 
Cartesian,  cylindrical,  spherical  or  generalized  nozzle  geometries,  depending 
on  the  value  of  a  logical  switch.  Since  the  algorithm  is  one-dimensional, 
timestep  splittina  in  the  various  directions  is  used  to  construct  two-  and 
three-dimensional  codes.  The  two-dimensional  Cartesian  and  cylindrical  cal¬ 
culations  are  actually  performed  with  the  same  code,  but  with  the  particular 
type  of  geometry  specified  at  the  beginninq  of  the  calculation. 

In  the  calculations  presented  below,  the  qrid  spacing  was  set  up  at  the 
beginning  of  each  calculation  and  held  fixed  in  time.  In  general,  the  cell 
spacing  should  not  change  more  than  20-30%  from  cell  to  cell.  For  the 
Cartesian  calculations  used  to  model  the  splitter  plate  experiments,  finely 
spaced  cells  were  clustered  around  the  centerline  where  the  instability  first 
occurs  and  the  coherent  structures  form.  For  the  cylindrical  calculations 
used  to  model  the  gas  jet,  the  grid  was  finely  spaced  in  the  jet  and  throuqh 
the  reqion  of  the  shear.  Sample  grids  are  shown  in  Figures  1  and  2. 

In  the  calculations  presented  below,  the  gravity,  diffusion,  chemistry, 
and  enerqy  release  options  included  in  the  code  are  not  used.  Thus  we  will 
not  describe  them  in  great  detail  in  this  report.  However,  we  note  that  the 
algorithms  for  chemistry,  energy  release  and  diffusion  have  been  discussed  by 
Oran  and  Boris  (1981). 
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IV.  OUTFLOW  BOUNDARY  CONDITIONS 


Solving  fluid  dynamics  problems  with  realistic  outflow  boundary  condi¬ 
tions  has  always  been  difficult.  The  fundamental  problem  is  that  information 
about  the  flow  beyond  the  computational  mesh  is  required  to  make  the  fluid 
near  the  boundaries  behave  properly.  There  are  a  number  of  ways  to  handle 
this  problem,  and  they  generally,  like  FCT,  involve  using  guard  cells  which 
are  not  actually  part  of  the  calculation  but  which  tell  the  boundary  cells  of 
the  computational  mesh  how  the  outside  world  is  behaving.  The  simplest  model 
of  outflow  in  guard  cells  is  to  say  that  the  momentum,  energy,  and  density  do 
not  change,  i.e.,  there  is  effectively  zero  gradient.  This  can  cause  prob¬ 
lems  in  long  time  calculations  since  it  does  not  provide  for  the  fact  that  as 
we  cro  further  from  the  phenomena  being  computed,  the  system  relaxes  to  back¬ 
ground  conditions. 

Here  we  present  an  outflow  algorithm  developed  to  use  with  the  FCT 
aloorithm  described  in  the  last  section.  The  outflow  algorithm  incorporates 
the  recruirement  that  the  solutions  must  relax  toward  ambient  conditions. 

Then  the  strong  nonlinear  stabilizing  properties  of  the  FCT  method  appear  to 
eliminate  instabilities  which  occur  in  other  nonlocal  methods  when  low  order 
extrapolations  are  used  for  specifying  boundary  conditions  (Turkel,  1980). 

The  types  of  fluid  problems  for  which  the  outflow  algorithm  is  intended  can 
be  written  as  a  coupled  series  of  single  or  multi -dimensional  continuity 
equations  of  the  form 


dg 

dt 


V»(Qv)  +  Source  -  Sinks 


(1) 
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where  the  compressible  fluid  must  somewhere  flow  off  the  edge  of  the  computa¬ 
tional  region.  In  this  equation  Q  is  a  conserved  quantity  such  as  the  mass. 


momentum,  or  energy  and  v  is  the  vector  fluid  velocity.  The  major  require¬ 
ment  of  the  algorithm  is  that  the  fluid  must  be  represented  in  the  region  off 
the  grid  by  fictitious  "guard"  or  "ghost"  cells.  These  values  are  used  to 
compute  the  derivatives  and  gradients  in  the  computational  cells  which  are  on 
the  edge  of  the  grid.  Throughout  the  discussion,  we  let  the  last  cell  on  the 
grid  be  cell  N.  The  subscript  g  will  be  used  to  indicate  the  guard  cell. 

The  linear  extrapolation  for  outflow, 

(2> 

is  unstable  when  used  in  conjunction  with  most  numerical  methods  for 
advancing  the  grid  variables.  Even  the  zeroth-order  extrapolation, 

?,*<>»  (3> 

is  unstable  in  the  sense  that  the  algorithm  has  no  knowledae  of  the  physi¬ 
cally  correct  asymptotic  value  of  Q,  and  hence  the  flow  cannot  ever  be 
expected  to  relax  to  ambient  conditions. 

The  formula  tested  here  is  a  simple  zeroth-order  extrapolation  plus  a 

slow  local  relaxation  toward  the  known  ambient  value  Q  .  . 

amb 


Q  ♦ 

VN  L 


char 


At 


char 


^  -  V 


(4) 


Here  L  .  is  a  characteristic  scale  for  the  flow  causing  the  relaxation  to 
char 

ambient,  and  vchar  is  the  velocity  of  this  relaxation.  The  Quantity  vchar 
is  typically  the  local  sound  speed  for  the  pressure  or  flow  velocity  for  an 
entropy  or  species  variable.  This  expression  is  an  approximation  to  the 
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lowest  order  terms  in  an  asymptotic  expansion.  It  is  valid  for  times  long 


compared  to  a  sonic  transit  time  of  the  system. 

Two  types  of  test  calculations  are  presented:  1)  a  diaphragm  breakina  in 
a  barrel  which  releases  a  high  pressure,  supersonic  gas  into  a  low  pressure 
background,  and  2)  a  cylindrical  jet  in  which  a  lower  density  gas  flows  into 
a  quiescent  background.  The  validity  of  the  algorithm  is  shown  by  comparing 
nested  series  of  calculations  where  the  outflow  boundary  of  one  grid  is 
interior  to  the  flowfield  of  a  larger  grid.  The  first  problem,  the  diaphragm 
in  a  barrel,  is  not  a  shear  flow  problem.  However,  it  is  included  in  this 
paper  because  it  is  an  initial  test  of  the  outflcw  algorithms  used. 


V.  TEST  OF  THE  OUTFLOW  ALGORITHM:  DIAPHRAGM  IN  A  BARREL 


1 .  Hie  Teat  Problem 

Below  we  compare  a  nested  series  of  calculations  where  the  outflow 
boundary  of  one  grid  is  interior  to  the  flew  field  of  a  larger  grid.  Then 
the  approximate  solution,  obtained  using  the  outflow  boundary  condition,  can 
be  calibrated  against  the  "correct"  solution  obtained  before  spurious 
information  from  the  outer  computational  boundary  arrives  at  the  outer  edge 
of  the  smaller  grid.  Calculations  with  three  different  size  grids  are 
presented,  a  40  x  40  grid,  an  80  x  80  grid,  and  an  even  larger  150  x  300 
grid. 

A  set  of  four  runs  using  the  FAST2D  model  have  been  performed.  The 
problem  chosen  is  cylindrically  symmetric  and  simulates  the  fluid  flows  which 
occur  when  a  thin  diaphragm  bursts.  The  diaphraom  initially  confines  a  high 
pressure,  isothermal  gas  in  a  thick-walled  barrel.  The  general  properties  of 
the  four  different  types  of  computations  performed  are  given  in  Table  I. 


Table  I. 


Calculation 

Grid  Size 

Boundary  Conditions 

Comments 

41 

40  x  40 

Reflecting  walls 

# 2 

40  x  40 

Outflow  Algorithm 

Test  of  new 

outflow  algorithm 

#3 

80  x  80 

Outflow  Algorithm 

Calibration 

for  Run 

#2 

#4 

150  x300 

Reflecting  walls 

Calibration 

for  Run 

#3 
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The  results  discussed  in  this  section  are  illustrated  by  digitized 
contour  plots  of  the  density.  Consider  Figure  3,  which  is  a  schematic 
diagram  of  the  initial  conditions  for  the  test  problem  considered  in  this 
section.  In  this  and  the  figures  following,  different  typed  letters  in  the 
computational  cells  combine  to  form  a  'shadowplot'  where  the  boundaries 
between  two  different  characters  represent  contours  of  the  fluid  density. 
Table  II  lists  the  letters  used  and  the  density  ranges  they  span  in  all  the 
figures  in  this  Section.  The  letter  "I"  is  used  for  alternate  bands  to 
improve  the  resolution. 

Table  II. 

Shadowplot  Contour  Levels:* 


1 .00E-03 

< 

< 

1  . 45E-03 

1 . 45E-03 

< 

- 

< 

2.10E-03 

2.10E-03 

< 

+• 

< 

3.20E-03 

3.20E-03 

< 

* 

< 

4.80E-03 

4.80E-03 

< 

A 

< 

7.00E-03 

7.00E-03 

< 

I 

< 

1  .  00E-02 

1 .00E-02 

< 

B 

< 

1  .45E-02 

1 .45E-02 

< 

I 

< 

2.10E-02 

2.10E-02 

< 

C 

< 

3.20E-02 

3.20E-02 

< 

I 

< 

4.80E-02 

4.80E-02 

< 

D 

< 

7.00E-02 

7.00E-02 

< 

I 

< 

1 .00E-01 

1 .00E-01 

< 

E 

< 

1 .45E-01 

1 .45E-01 

< 

I 

< 

2.1 OE-01 

2.1 0E-01 

< 

p 

< 

3.20E-01 

3.20E-01 

< 

I 

< 

4.80E-01 

4.80E-01 

< 

G 

< 

7.00E-01 

7.00E-01 

< 

I 

< 

1 .00E+00 

1 .00F+00 

< 

H 

< 

1 .45E+00 

1 .45E+00 

< 

I 

< 

2.10E+00 

read  "E-03" 

as 

10' 

-3 
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The  initial  conditions  are  the  same  for  calculations  #1-4.  The  ambient 
density  is  1.29  *  10"3  g/cm"3,  the  pressure  is  a  fixed  multiple,  1.0  x  io9, 
of  the  density,  and  the  lower  half  of  the  barrel  is  pressurized  to  1000  times 
the  ambient  pressure.  The  barrel  geometry  i3  the  same  for  all  cases  listed 
in  Table  I.  The  spacing  of  the  computational  grid,  0.1  cm,  is  also  fixed  in 
all  the  computations.  The  cases  listed  in  Table  I  differ  only  in  the  loca¬ 
tion  of  the  upper  and  right  hand  boundary  and  in  the  boundary  condition 
applied . 


2.  The  Outflow  Algorithm 

The  conserved  quantities  were  convected  by  FCT  in  this  problem:  the 
mass  density,  p,  and  the  two  components  of  the  momentum  density,  pvx  and 
pvy.  The  characteristic  velocity  for  relaxation  of  density  and  pressure 
fluctuations  is  the  sound  speed,  here  constant  at  3.2  x  io4  cm/s.  The  char¬ 
acteristic  lenqth  was  taken  as  the  radius  of  the  cylinder,  4.0,  3.0,  and  15.0 
cm  for  the  three  different  sized  grids.  The  density  was  extrapolated  beyond 

the  end  of  each  exterior  row  or  column  by 

C 


pN*ir  * 

wall 


(Pamb  -  V 


(5) 


(pVg  3  (pVn  (6 

{pvz)g  "  (pv,N  (7 

where  Cg  is  the  sound  speed.  Since  the  pressure  and  density  are 

proportional  to  each  other,  the  pressure  is  extrapolated  as  in  Eg.  (5). 
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3.  Results  of  the  Calculations 


Figure  4  shows  a  series  of  shadowplots  from  Calculation  #1 ,  which  has 
reflecting  walls.  Figure  4A  shows  that  after  25  timestepe  the  shock  has 
moved  only  about  5  cells  up  the  barrel  (0.5  cm)  and  has  not  yet  reached  the 
rim.  By  step  50,  shown  in  Figure  4B,  the  flow  inside  the  barrel  is  still 
essentially  one-dimensional.  However,  as  the  shock  and  following  supersonic 
flow  emerge  from  the  barrel,  the  first  evidence  of  the  multidimensional 
nature  of  the  flow  is  evident.  By  step  100,  shown  in  Figure  4C,  the  flow  has 
hit  the  upper  boundary  and  the  rarefaction  is  working  its  way  back  toward  the 
bottom  of  the  barrel.  By  step  200,  shown  in  Figure  4D,  the  primary  shock  has 
had  time  to  rattle  bade  and  forth  between  the  solid  barrel  wall  and  the  outer 
and  top  boundaries,  which  here  are  solid  walls. 

Figure  5  presents  a  series  of  shadowplots  taken  from  Calculation  #4, 
performed  on  the  150  x  300  grid.  These  plots  show  the  solution  of  the  out¬ 
flow  test  problem  without  any  interference  from  the  boundaries.  The  barrel 
wall  and  the  nested  40  x  40  and  80  *  80  computational  regions  are  outlined  on 
this  figure.  Note  that  of  the  150  x  300  cells  used  in  the  computation,  only 
128  x  215  are  shown  on  the  figures. 

Figures  6  through  10  are  composites  of  the  different  calculations  at 
fixed  timesteps.  Calculations  *2  and  3  have  the  outflow  boundary  conditions 
at  the  upper  and  right  hand  walls.  From  Figure  6  -we  3ee  that  by  step  100  the 
shock  would  have  reflected  from  the  top  wall  in  Calculation  #1.  Comparing 
Figures  4  and  6,  we  see  that  by  this  time  the  reflection  from  the  closed  wall 
in  Figure  4  has  affected  the  answers  down  to  about  cell  34  on  the  Z  axis. 
Within  and  just  outside  of  the  barrel,  the  solutions  are  still  essentially 
the  same.  Calculation  #3  in  Figure  6  shows  that  the  expanding  plume  has 


reached  cell  65  (Z  ■  6.5  cm)  and  so  the  flow  is  still  interior  to  the  80  x  80 


grid.  At  this  point  no  errors  have  propagated  back  from  this  boundary  into 
the  interior  region.  The  outline  of  the  region  considered  in  Calculations  #1 
and  #2  are  shown  in  this  figure.  Calculation  #3,  which  uses  the  80  x  80 
grid,  and  the  new  outflow  boundary  conditions,  provides  a  calibration  of  the 
open  boundary  condition  used  in  Calculation  #2.  Thus  we  can  here  compare 
Calculations  42  and  43  in  Figure  6  and  see  that  they  show  substantially  the 
same  result. 

In  Figure  7  (step  150),  Calculation  #3,  the  rarefaction  behind  the 
primary  shock  is  well  formed  and  appears  as  a  white  strip  starting  from  the 
outer  corner  of  the  barrel  wall.  Reflection  off  the  upper  boundary  at  Z  = 
8.0  cm  would  have  occurred  by  this  time  if  the  upper  wall  had  been  a  reflec¬ 
ting  boundary.  In  this  case,  where  we  are  using  a  model  for  the  outflow,  we 
can  see  that  small  errors  are  beginning  to  move  inward. 

By  step  200  in  Figure  8,  Calculation  #3,  the  flow  has  exited  through  the 
upper  surface  but  has  not  yet  reached  the  outer  wall.  Some  errors  have  been 
propagating  toward  the  region  of  the  40  x  40  calculation  on  the  larger 
80  x  80  grid,  but  have  not  reached  it  yet.  Thus  this  Figure  is  good  for 
comparison  with  Figure  8,  Calculation  #2.  The  extent  of  the  outflow  boundary 
condition  errors  at  the  boundary  of  the  80  x  80  region  is  found  by  comparing 
with  Figure  8,  Calculation  #4. 

By  step  300  in  Figure  9,  Calculation  #2,  the  primary  shock  has  spread 
down  around  the  barrel  and  has  reached  the  ground  plane  at  Z  »  0.0  cm.  The 
errors  which  are  generated  by  the  open  boundary  conditions  should  have 
reached  the  40  x  40  region  by  now,  yet  Calculations  #2,  3,  and  4  are  still 
essentially  identical. 
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Finally,  by  step  400  in  Figure  10,  Calculation  #3,  the  errors  from  the 


outflow  boundary  conditions  have  reached  the  vicinity  of  the  barrel.  Compar¬ 
ison  with  Figures  10,  Calculations  #2  and  4,  show  only  small  differences 
which  can  be  associated  with  the  outflow  boundary  algorithm  we  are  testing. 
The  accuracy  of  this  result  in  the  40  *  40  case  is  surprising  where  the 
computational  boundary  was  closer  than  a  barrel  diameter  at  every  point. 
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VI.  TEST  OF  THE  OUTFLOW  ALGORITHM:  CYLINDRICAL  GAS  JET 


1 .  The  Test  Problem 

The  test  problem  described  above  is  related  to  interior  ballistics, 
atmospheric  explosions  and  airblast  problems  as  well  as  shear  flows.  It 
shows  our  first  application  of  these  boundary  conditions  in  an  important 
problem  where  both  supersonic  and  subsonic  flows  are  important.  The  flows 
were  dominated  by  divergence  rather  than  rotation,  and  the  relaxation  was 
primarily  by  acoustic  waves  rather  than  convection  of  vortices  off  the 
system  boundary.  The  excellent  results  obtained  from  these  tests  encouraged 
us  to  proceed  with  the  difficult  subsonic  problems  described  here  in  this 
Section.  In  combustion  problems,  expansion  and  vorticity  are  equally 
important.  Thus  we  felt  that  the  separate  tests  of  both  fluid  dynamic 
aspects  was  appropriate. 

The  Kelvin-Helmholtz  outflow  tests  presented  in  this  section  were 
done  in  preparation  for  transition  to  turbulence  calculations  such  as  those 
presented  in  the  next  section.  We  wish  to  test  the  boundary  condition 
algorithm  described  above  in  a  calculation  in  which  the  flows  are  subsonic 
throughout  the  course  of  the  calculation  and  hence  the  flow  is  essentially 
incompressible.  Here  we  consider  a  jet  of  neon  gas  moving  at  2  *  104  cm/s 
which  is  vented  into  a  quiescent  air  background.  There  is  a  small  initial 
perturbation  at  the  interface  between  the  jet  and  the  air  which  causes  the 
the  shear  flow  eventually  to  become  unstable.  After  a  longer  time,  the  usual 
pattern  of  rotating  coherent  structures  forms  at  the  interface.  These 
structures  merge  and  flow  out  of  the  computational  system,  i.e.,  off  the  top 
of  the  computational  mesh. 
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2.  The  Outflow  Algorithm 

The  basic  outflow  algorithm  used  in  this  case  states  that  the  pressure 

in  the  guard  cell  just  across  the  boundary  is  extrapolated  from  the  boundary 

cell  but  relaxes  to  a  prespecified  ambient  value.  Specifically, 

P  =  P„  +  t  •  (P  .  -P„)  (8) 

a  N  amD  N 


where  P  is  the  guard  cell  pressure,  p  ^  is  the  ambient  pressure,  and  P .  is 
g  ‘  a mb  N 

the  boundary  pressure.  Here  ui  reflects  the  rate  of  relaxation  and  is  con¬ 
trolled  by  the  sound  speed,  Cg,  and  a  characteristic  radius,  R* .  We  have 
taken  R'  as  the  characteristic  radius  of  the  disturbance,  in  this  case  the 
radius  of  the  nozzle.  Since  the  flows  here  are  typically  a  third  to  a  fourth 
of  the  speed  of  sound,  the  outflow  velocity  for  vorticity  can  be  estimated 
using  the  sound  speed.  For  slower  flows  the  fluid  velocity  should  be  used. 
The  guard  cell  pressure.  Equation  (85,  then  effects  the  calculation  of  the 
guard  cell  energy  value.  The  relaxation  rate  a»  is  an  adjustable  parameter 
representing  the  characteristic  rate  and  size  of  flew  patterns  leaving  the 
grid  and  presumably  could  be  adjusted  to  improve  the  match  between  large-  and 
small-grid  calculations  still  further.  The  density  and  momenta  are  extrapo¬ 
lated  according  to 


(  pv  )  =  (  pv  ) 

x  a  x 


(PVg  +  (pyN 
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Note  atht  in  this  case  unlike  the  barrel  caclulation  previously  described, 
pressure  and  density  are  decoupled.  Thus  to  account  for  effects  of  sound 
waves,  we  need  to  modify  only  the  pressure  at  the  boundary.  Modifying  the 
density  has  little  or  no  effect  because  there  is  no  substantial  inflow  of 
mass  except  at  the  nozzle. 

3.  Results  of  the  Calculations 

In  order  to  test  the  outflow  boundary  conditions  algorithm,  we  performed 
a  series  of  three  calculations  summarized  in  Table  III.  In  the  first  and 
largest  calculation,  a  stretched  60  *  120  mesh  represents  a  region  7.6  cm  in 
radius  by  26.0  cm  in  length.  The  second  calculation  was  initialized  by  using 
the  results  from  step  2000  of  the  first  calculation,  but  uses  only  the  first 
54  x  100  cells  of  the  original  mesh.  Thus  it  represents  a  volume  of  3.25  cm 
radius  by  10.8  cm  length.  The  final  calculation  used  a  54  *  100  mesh, 
equivalent  to  the  mesh  used  in  Calculation  #2.  However,  this  last 
calculation  was  started  from  step  =  0  with  the  smaller  mesh. 

Table  III 


Calculation 

Grid  Size 

Comments 

41 

60  x  120 

Step  0  -  3000 

42 

54  x  ioo 

Steps  2000  -  3000 

Restarted  from  step  2000  of  #1 

#3 

54  x  ioo 

Steps  0  -  3000 
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In  order  to  compare  the  results  of  the  two  calculations,  we  present 


contours  of  the  quantity 


R  = 


n  ( 1  ) 


(13) 


n(  1 )  +  n(  2) 

where  n( 1 )  is  the  number  density  of  neon  and  n(2)  is  the  number  density  of 
air,  where  air  was  taken  as  one  species  with  the  mean  properties  of  a  mixture 
of  20%  02  and  80%  N2 .  By  definition,  the  contours  range  from  R  =  0,  100% 
air,  to  R  =  1,  100%  neon.  Figure  11  shows  the  ratio  R  for  Calculation  #1  at 
steps  500,  1000,  2000,  and  3000  and  presents  a  time  history  of  the  flow.  We 
see  that  the  system  initially  goes  unstable  at  1.14  cm  from  the  nozzle.  The 
structures  grow,  merge  and  finally  exit  out  of  the  top  of  the  computational 
grid.  We  have  set  the  inflow  velocity  to  2  *  10 4  cm/s.  It  takes  ~1000  time 
steps  for  a  fluid  element  entering  the  system  at  the  nozzle  to  reach  the  top 
of  the  large  arid,  26  cm  if  it  is  not  slowed  appreciably  in  interaction  with 
the  background  air.  Thus  much  of  the  material  entering  at  step  2000  will  be 
exiting  the  computational  arid  shortly  after  step  3000.  The  material  effec¬ 
tively  chances  over  every  1000  steps.  Thus  the  calculation  was  continued  far 
enouch  to  test  the  boundary  condition  algorithm  for  long  term  fidelity  and 
stability. 

Figure  12  shows  a  comparison  of  step  3000  for  calculations  #1,  2,  and  3. 
First  we  note  that  the  first  two  figures  look  very  similar  within  the 
54  x  100  range,  which  is  very  encouragina.  The  last  figure,  on  the  far 
right,  however,  while  qualitatively  similar  to  the  first  two,  has  differ¬ 
ences.  Figure  13  shows  a  composite  for  the  54  x  100  and  60  *  120  calcula¬ 
tions,  starting  at  step  500  through  step  1500.  Through  steps  1000  the  calcu¬ 
lations  are  quite  similar,  however  small  but  noticeable  differences  begin  to 
appear  between  steps  1000  and  1500. 
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4.  Discussion 


The  constants  chosen  in  implementing  the  boundary  conditions  may  not 
have  been  the  best  ones  and  this  has  to  be  tested.  First,  the  value  of  R' 
should  probably  be  larger  than  the  radius  of  the  inlet.  This  value  should 
reflect  the  characteristic  sizes  of  the  structures,  and  we  knew  that  they 
grow  according  to  a  similarity  condition.  Also,  the  appropriate  velocity  to 
use  would  be  the  fluid  velocity,  not  the  speed  of  sound,  if  the  vortices  were 
leaving  the  region  at  very  slow  speed.  The  speed  of  sound  was  appropriate  in 
the  barrel  problem  because  that  was  a  problem  in  supersonic  flow.  Here  we 
are  concerned  with  subsonic  flows. 

The  shear  flow  discussed  in  this  section  is  unstable  and  will  always 
go  unstable  as  long  as  there  is  some  small  perturbation,  even  just  noise  or 
roundoff,  between  the  flowing  and  quiescent  gas.  However,  unlike  the  dia- 
phracim  in  the  barrel  problem  discussed  above,  the  exact  appearance  of  the 
Kelvin-  Helmholtz  instability  in  the  nonlinear  regime  is  very  sensitive  to 
background  conditions,  initial  conditions,  and  to  small  f luctuat.  ons  in  the 
system.  This  is  the  crux  of  the  problem  in  experiments  a\Ro:  many  realiza¬ 
tions  are  possible  depending  on  fluctuations  in  the  initial  and  boundary 
conditions.  Tn  these  situations  some  of  the  difficulty  encountered  in  doing 
the  computations  reflects  the  same  difficulty  in  the  experiments.  Both  high¬ 
light  the  highly  nonlinear  sensitivity  of  the  system  to  the  boundary  condi¬ 
tions  . 
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VII.  GAS  JET  AND  SPLITTER  PLATE  SIMULATIONS 

In  this  section  we  describe  preliminary  results  obtained  using  the  model 
described  above  to  simulate  the  coherent  structures  developed  from  Kelvin- 
Helmholtz  instabilities  in  shear  flows.  In  the  previous  section,  we  empha¬ 
sized  numerical  methods;  here  we  emphasize  the  physical  results.  We  discuss 
two  configurations  of  the  FAST2D  model  discussed  above:  a  cylindrical  gas  jet 
and  the  Cartesian  splitter  plate  geometry. 

1 .  Splitter  Plate 

Figure  14  shows  a  sequence  of  frames  from  a  two -dimensional  simulation 
of  the  flew  generated  at  an  idealized  splitter  plate.  The  tip  of  the  split¬ 
ter  plate  is  located  at  x  *  0.0  cm  and  y  =  4.0  cm,  the  dividing  point  between 
the  light  gray  fluid  entering  into  the  lower  half  of  the  computational  region 
(air  at  high  velocity)  and  the  darker  gray  fluid  (air  at  somewhat  lower 
velocity) .  Figure  2  shows  the  type  of  variably  spaced  grid  used  near  the 
trailing  edge  to  represent  the  Cartesian  splitter  plate  problem.  The  y- 
direction  has  good  resolution  near  the  centerline.  The  velocity  of  the 
lower  stream  is  1.5  x  104  cm/s  and  the  upper  stream  is  considerably  slower, 
flowing  from  left  to  right  at  5  x  10 3  cm/s.  As  the  large  spanwise  vortices 
propagate  downstream,  they  merge  and  grow.  Therefore  the  larger  computa¬ 
tional  cells  encountered  by  the  flow  further  downstream  and  further  from  the 
centerline  will  be  adequate  to  resolve  the  larger  coherent  structures  about 
as  well  as  the  smaller  structures  are  resolved  near  the  splitter  plate. 

The  initial  stages  of  the  development  of  the  flow  are  not  shown  in 
Figure  14.  During  this  early  period,  before  0.87  ms,  the  shear  flow  estab¬ 
lished  by  the  trailing  edge  of  the  splitter  plate  is  initially  restricted  to 
a  very  narrow  layer  which  broadens  quickly  at  first  due  to  viscosity  and  then 
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more  slowly  as  both  streams  of  fluid  leave  the  splitter  plate.  As  in  the 
case  of  the  cylindrically  symmetric  jet  simulations  discussed  in  the  boundary 
condition  test  problem  and  in  the  next  section  below,  the  flow  is  Kelvin- 
Helmholtz  unstable.  Initially,  during  the  linear  growth  period  of  the  insta¬ 
bility,  most  of  the  computed  region  behaves  as  if  there  were  periodic  bound¬ 
ary  conditions  because  the  impressed  perturbation  has  a  wavelength  short 
compared  to  the  system  length.  By  0.87  ms,  the  first  panel  in  Figure  14,  the 
influence  of  the  inflow  and  outflow  is  felt.  Merging  vortices  move  off  the 
right  edge  of  the  computational  system  and  unperturbed  fluid  moves  in  on  the 
left. 

Our  preliminary  calculations  used  boundary  conditions  which  prespecified 
values  for  the  mass,  momentum  and  energy  density  of  the  inflowino  gas.  These 
conditions,  however,  did  not  correctly  provide  the  feedback  from  the  Kelvin- 
Helmholtz  rolls  and  vortex  mercting  on  to  the  newly  entering  fluid.  We  know 
that  pressure  pulses  from  upstream  create  small  transverse  flows  at  the 
trailing  edge  of  the  splitter  plate.  These  pulses,  then,  start  the  insta¬ 
bility  at  finite  amplitude  for  the  next  coherent  vortex  roll  up.  when  the 
inflow  pressure  was  specified  as  ambient,  we  observed  that  the  first  vortex 
structures  formed  very  far  downstream.  The  apparent  reason  for  this  was  that 
the  pressure  perturbations  arriving  at  the  inflow  boundary  were  cancelled  by 
the  non-physical  condition  of  specifying  the  inflow  pressure.  Clearly  a  more 
physically  reasonable  treatment  of  the  inflow  response  to  pressure  fluctua¬ 
tions  as  well  as  the  outflow  was  required  for  this  type  of  problem. 

The  more  correct  way  to  treat  the  inflow  is  to  specify  the  inflow 
density  and  velocity,  and  then  to  use  a  zero  slope  condition  on  the  pressure 
at  the  inflow  boundary  to  derive  the  energy.  This  algorithm  was  used  for  all 
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the  simultions  which  follow,  it  allows  pressure  differences  between  the  top 
and  bottom  streams  to  generate  transverse  flows  and  implies  that  the  physical 
plenum  admitting  the  gas  into  the  region  we  are  simulating  is  of  zero  length. 
Physical  treatment  of  a  real  plenum  could  be  added  by  extending  the  splitter 
plate  into  the  computational  region.  A  simple  lumped  parameter  boundary 
calculation  representing  the  real  plenum  could  also  be  developed  using  the 
more  detailed  model  to  calibrate  it.  Then  a  phenomenology  would  be  incorpor¬ 
ated  to  estimate  the  pressure  build  up  and  velocity  changes  which  could  be 
expected  on  inflow.  Thus  we  see  that  the  pressure  reflection  condition  is 
one  limit  of  such  a  model  for  which  the  plenum  has  zero  volume.  This  limit 
is  appropriate  when  the  fluid  flow  is  slow  compared  to  the  sound  speed,  the 
situation  most  often  considered  in  experiments  but  not  necessarily  applicable 
when  the  flows  are  as  fast  as  in  this  calculation. 

The  outflow  boundary  conditions  were  handled  as  described  in  the  previ¬ 
ous  sections  of  this  paper.  With  a  fluctuating  inflow  pressure  it  is  import¬ 
ant  to  relax  the  outflow  pressure  toward  an  ambient  value.  This  is  because  a 
background,  base  pressure  for  the  system  has  to  specified  in  compressible 
calculations  when  the  inflow  pressure  is  allowed  to  float  in  response  to 
pressure  waves  from  upstream.  Relaxation  of  the  density  and  momentum  toward 
ambient  is  less  important  on  outflow  because  these  values  are  given  at  the 
inflow  boundary.  The  long  time  values  of  all  the  variables  must  be  available 
to  the  solution  from  the  beoinning  to  prevent  secular  deviation  of  the  calcu¬ 
lated  flow  from  properties  characterizing  the  ambient  medium. 

The  mechanism  which  reinitiates  the  instabilities  close  to  the  inflow 
boundary  works  through  pressure  pulses  generated  at  various  scales  by  the 
fluid  accelerations  involved  in  nonlinear  vortex  rollup  and  vortex  merging 
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downstream.  They  are  required  to  ensure  that  the  flow  remains  essentially 
diverqence  free  everywhere  when  a  coherent  structure  or  vortex  is  locally 
accelerated.  Though  these  pulses  are  transmitted  acoustically,  they  exist 
even  in  the  incompressible  limit. 

Evidence  for  this  mechanism  is  supplied  by  the  sequence  of  density 
shadowplots  shown  m  Ficrure  14.  By  0.87  ms  the  initial  linear  growth  phase 
is  over  and  the  instability  has  saturated  rather  uniformly  along  the  original 
slip  line.  States  similar  to  what  occurred  in  this  calculation  prior  to  0.87 
ms  are  shown  in  the  following  section  on  the  evolution  of  the  instability  in 
the  round  jet.  In  the  splitter  plate  calculation,  a  long  unrippled  portion 
of  the  interface  is  seen  on  the  left  of  the  0.87  ms  panel,  showing  that  the 
initial  instability  has  flowed  away  from  the  trailing  edge.  However,  the 
pressure  perturbations  resulting  from  this  were  small  and  short  scale,  so 
that  they  did  not  initiate  the  Kelvin-Helmholtz  instability  strongly  enough 
for  small  scale  vortices  to  appear  on  the  left  of  the  0.87  ms  panel. 

The  first  vortex  merging  is  also  shown  approaching  the  rioht  hand 
boundary  in  the  top  three  panels  of  the  figure.  This  flow  structure  is  both 
large  and  irregular  enough  to  Generate  stronger  pressure  perturbations,  which 
can  be  seen  driving  noticeable  fluctuations  into  the  left  half  of  the  inter¬ 
face  at  0.98  and  1  .03  ms.  By  the  last  frame  at  2.06  ms,  a  spectrum  of  modes 
is  seen  in  what  looks  like  a  snapshot  from  an  experimental  flow  visualiza¬ 
tion. 

2.  Gas  Jet 

First  we  consider  the  coherent  structures  in  the  transition  from  laminar 
to  turbulent  flow  in  coflowing  gas  jets  in  two  cases:  air  into  air,  which  is 
constant  density,  and  air  into  freon,  where  the  density  ratio  is  about  6:1. 
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These  simulations  were  performed  in  order  to  isolate  and  study  the  effects 
of  density  gradients,  which  is  one  of  the  the  major  physical  conditions 
distinguishing  reactive  flow  from  classical  turbulence.  The  calculations 
also  provided  important  information  on  the  application  of  inflow  and  outflow 
boundary  conditions  to  subsonic  shear  flow  problems. 

In  the  first  case,  air  into  air,  the  densities  are  constant  and  these 
density  gradients  do  not  exist.  However,  the  second  case,  air  into  freon,  a 
light  gas  into  a  heavy  gas,  has  the  added  effect  of  vorticity  generation 
through  the  term 

„4, 


where  5  is  the  vorticity,  and  o  and  P  are  the  fluid  density  and  pressure 
fields,  respectively. 

The  geometry  is  the  same  as  in  the  gas  jet  described  in  Section  VI.  The 
air  flows  into  the  cylinder  from  a  hole  at  the  bottom  center  at  a  velocity  of 
1.5  x  104  cm/s.  Unlike  the  problem  discussed  in  Section  VI,  in  which  the  jet 
flowed  into  a  quiescent  background,  the  fast  jet  of  air  is  surrounded  by  a 
co-flowing  jet  at  the  much  lower  velocity  of  5  x  10 3  cm/s.  The  co-flow  is 
either  air  or  freon.  Also,  the  right  hand  boundary  here  is  a  hard  wall  with 
only  the  top  boundary  open. 

The  instability  in  the  fluid  is  generated  by  imposing  a  small  sinusoidal 
perturbation  on  the  momentum  at  the  interface,  initially  at  r  *  0.6  cm.  This 
divergence  free  perturbation  is  largest  at  the  material  discontinuity  and 
goes  smoothly  to  zero  at  the  sidewalls  and  at  the  center  of  the  system.  The 
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maximum  amplitude  of  the  perturbation  is  2%  of  the  jet  velocity.  Given  the 
initial  density  and  the  momentum  perturbation,  an  incompressible  velocity 
perturbation  is  implied. 

Figures  15  and  16  summarize  the  air  into  air  and  the  air  into  freon 
calculations  we  wish  to  discuss.  The  first  effect  observed  in  these  calcula¬ 
tions  is  the  formation  of  a  Kelvin-Helmholtz  instability  at  the  3lip  line. 
Because  of  the  short  wavelength  and  fixed  amplitude  of  the  perturbation,  the 
problem  initially  evolves  as  a  periodic,  time  dependent  system  except  for  one 
or  two  instability  wavelengths  at  each  end.  The  first  panel  of  Figure  15 
(the  air  into  air  cylindrical  jet)  shows  this  effect  at  240  us.  Here  a 
rather  uniform  band  of  mixed  material  forms  separating  the  two  fluids.  Thi3 
uniform  buffer  zone  does  not  initially  display  the  vortex  meroing  which  char¬ 
acterizes  the  fully  developed  nonlinear  flow.  Indeed,  a  fully  developed 
turbulent-looking  flow  is  not  set  up  until  most  of  the  jet  material  initially 
on  the  computational  mesh  has  flowed  out  the  top  and  been  replaced  by  fresh 
material  through  inflow. 

'rtie  two  calculations  look  qualitatively  different.  There  is  a  notice¬ 
able  decrease  in  the  entrainment  in  the  air-into-f reon  calculation.  The 
vorticity  generated  at  the  nozzle  moves  away  at  a  density-weighted  velocity 
intermediate  between  the  velocity  of  the  two  streams.  In  the  air  into  freon 
case,  Figure  16,  the  sloped  lines  on  the  figure  which  track  structures  moving 
on  interface  shor  that  these  structures  move  more  slowly  than  in  the  air  into 
air  case  in  Figure  15.  As  the  coherent  structures  move  with  a  density 
weighted  velocity,  the  air  into  freon  case  has  a  much  greater  slip  velocity 
between  the  jet  and  the  coflow.  As  a  result,  for  the  air  into  freon  case, 
the  freon  entrainment  is  reduced  and  the  nonlinear  coherent  structures  form 
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much  closer  to  the  nozzle  lip.  As  can  be  seen  from  the  sloped  lines  on  the 
two  figures,  the  velocity  of  the  vortices  is  faster  in  the  air  into  air  case 
by  about  a  factor  of  two.  Because  the  air  is  moving  much  faster  than  the 
vortices  in  the  air  into  freon  case,  relatively  little  freon  gets  entrained. 

The  calculation  shown  in  Figure  16  shows  a  choking  effect  in  which  the 
dense  coherent  structures  squeeze  light  material  in  the  jet  and  push  the  jet 
flow  outward  at  the  nozzle.  This  aoain  affects  the  lower  boundary  through 
the  pressure  reflection  condition,  and  thus  effects  the  inflow  of  material. 
Then  the  Bernoulli  effect  accelerates  the  jet  fluid  through  the  slowly 
changing  channel  formed  by  the  slower,  higher  density  coherent  structures. 

The  panel  at  390  us  on  the  lower  left  shows  the  jet  accelerating  so  much  that 
the  density  has  dropped  noticeably.  A  higher  pressure  region  forms  just 
above  the  nozzle  to  accelerate  the  jet  axially  through  the  constriction 
formed  by  the  coherent  structure.  This  higher  pressure  also  pushed  the 
interface  laterally  at  the  nozzle  lip,  as  can  be  seen  in  the  remaining  panels 
of  Figure  16.  This  displacement  becomes  the  seed  for  the  next  coherent 
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VIII.  DISCUSSION  AND  SUMMARY 


this  paper  has  presented  a  summary  of  work  to  date  at  NRL  on  developing 
and  testing  a  model  to  use  for  studies  of  unstable  shear  flows  transitioning 
to  turbulence.  The  main  features  of  the  model  are  that  it  is  two-dimen¬ 
sional,  fully  compressible,  and  time  dependent  in  either  Cartesian  or  cylin¬ 
drical  geometry.  To  date,  emphasis  has  been  on  development  of  the  model  as 
opposed  to  application,  although  some  results  and  model  predictions  have  been 
discussed  in  Section  VII . 

Future  emphasis  will  be  on  applications  althouoh  there  is  is  still  much 
to  be  learned  about  inf lew  perturbations.  In  particular,  we  are  interested 
in  how  well  this  model  compares  to  experiments.  Calculations  using  this 
model  will  be  compared  to  experiments  to  test  the  effects  of  density  differ¬ 
ences  on  the  initial  instability  and  on  transition  to  turbulence.  We  now 
know  that  it  be  necessary  to  add  a  number  of  other  physical  processes  to  the 
purely  convective  calculations  shown  above.  For  example,  we  might  have  to 
consider  buoyancy  and  molecular  diffusion.  Buoyancy  should  be  important  at 
distances  8  to  10  nozzle  diameters  downstream,  and  molecular  diffusion  could 
be  important  when  very  light  gas  flows  into  a  heavier  gas.  Algorithms  for 
gravity,  and  molecular  diffusion  are  currently  a  part  of  the  model.  We  are 
currently  developing  subgrid  turbulence  models,  and  the  particular  inflow 
nozzle  model  will  have  to  be  developed  in  coordination  with  the  experiments. 
Finally,  computational  diagnostics  must  be  included  that  operate  on  the  com¬ 
puted  primitive  quantities  (such  as  mass,  momentum  and  energy)  and  put  them 
in  a  form  that  can  be  directly  compared  with  the  experiments.  These  diagnos¬ 
tics  includes  such  quantities  such  as  fluctuation  velocities,  correlations, 
and  various  averages. 
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We  will  also  look  at  the  effects  of  energy  release,  which  occurs  in 


reacting  shear  flows  in  combustors.  Calculations  with  a  model  for  heat 
release  should  at  least  qualitatively  reproduce  some  of  the  bulk  effects  seen 
in  the  experiments  (Yule  and  Chigier,  1979;  Yule,  1981).  For  example,  there 
is  a  noticeable  extension  of  the  turbulent  transition  region  when  there  is 
chemical  energy  release. 

In  conclusion,  we  feel  that  we  now  have  a  tool  which  has  been  developed 
carefully  and  tested  in  a  number  of  calculations.  Hie  limitations  and 
remaining  difficulties  are  reasonably  well  known. 
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Figure  2.  Sample  computational  grid  used  in  the  Cartesian  calculations  of 
the  splitter  plate. 
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Fiqure  3.  Initial  conditions  for  the  calculations  of  an  explosion  in  a 

barrel.  The  outer  and  upper  walls  are  solid  in  Calculations  #1 
and  #4.  Calculations  #2  and  3  use  the  outflow  algorithm  at  the 
upper  and  outer  walls. 
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Figure  4.  3arrel  explosion  problem  for  Calculation  #1  at  four  times  in  the 
calculation.  Reflection  from  the  top  and  right  hand  side  is 
observed  in  C  and  D  due  to  the  solid  walls. 


Fiyure  7.  Step  150  for  barrel  Calculations  #2 
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CALCULATION  *4 


Figure  10.  Step  400  for  barrel  era lculations  #2,  3,  and 


aumjiimiuLuiSooM  . . . . mi . .  Biimumniuiuinuuuu 


Comparison  of  contour  plots  of  the  mixing  parameter  F  at  step  300 
for  three  calculations:  60  *  120,  54  x  100  restarted  from  the 
large  calculation  at  step  2000,  and  a  54  *  100  calculation. 
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Fiqure  15*  Shadow  plots  of  the  mixing  pa 
fast  air  jet  moving  into  a  sli 


47 


m 


•  '  tJ 

;«  ‘j 

;  .  J 

1.  ■:S'S;' :;’  . 

•'.V-fW-i.  ‘ .  "  > 

REFERENCES 


Book,  D.,  J.  Boris,  A.  Kuhl,  E.  Oran,  M.  Picone,  S.  Zalesak  (1980), 

Simulation  of  Complex  Shock  Reflections  from  Wedge  in  Inert  and  Reactive 
Gas  Mixtures,  Proceedings  of  the  7th  International  Conference  on 
Numerical  Methods  in  FLuid  Dynamics,  pp.  84-90,  Sprincer-Verlag,  New 
York. 

Boris,  J.P.  (1976a),  Flux  Corrected  Transport  Modules  for  Solving  Generalized 
Continuity  Equations,  Naval  Research  Laboratory  Memorandum  Report  3237. 

Boris,  J.P.,  (1976b)  Numerical  Solution  of  the  Continuity  Equation, 

Proceedings  of  the  Second  European  Conference  on  Computational  Physics, 
North  Holland  Pub.  Co.,  1976.  Also,  NRL  Memorandum  Report  3327,  Naval 
Research  Laboratory,  Wash.,  D.C.,  20375. 

Boris,  J.P.,  and  D.L.  Bock  (1976),  Solution  of  Continuity  Equations  by  the 

Method  of  Flux  Corrected  Transport,  in  Methods  in  Computational  Physics , 
Vol .  16,  85-129,  Academic  Press,  New  York. 

Browand,  F.K.,  and  P.D.  Weidman  (1976),  Large  Scales  in  the  Developing 
Mixing-Layer,  J.  Fluid  Mech,  76,  127-144. 

Browand,  F.K.,  and  C.D.  winant  1973),  Laboratory  Observations  of  Shear  Layer 
Instability  in  a  Stratified  Fluid,  Boundary-Layer  Meteor.,  5,  67. 

Brown,  G.,  and  A.  Roshko  (1974),  On  Density  Effects  and  Larne  Structure  in 
Turbulent  Mixing  Layers,  J.  Fluid  Mech.  64,  775-816. 

Chigier,  N.A.,  and  A.J.  Yule  (1979),  The  Physical  Structure  of  Turbulent 
Flames,  AIAA  Paper  79-0217,  New  York. 

Landau,  L.D.,  and  E.M  Lifshitz  (1959),  Fluid  Mechanics,  Peraamon  Press,  New 
York .  ~~ 

Liepmann,  H.W.  (1979),  The  Rise  and  Fall  of  Ideas  in  Turbulence,  American 
-Scientist,  67,  221-228. 

Oran,  S.S.,  and  J.P.  Boris,  Detailed  Modellina  of  Combustion  System,  prog. 
Energy  Combustion  Sci .  7,  1,  1981. 

Oran,  E.S.,  T.R .  Young,  J.P.  Boris,  J.M.  Picone,  D.H.  Edwards  (1982) 

Nineteenth  Symposium  (International)  on  Combustion,  1982,  p.  573.  The 
Combustion  Institute,  Pittsburgh. 

Overman,  E.A.,  and  N.J.  Zabusky  (1981),  Transition,  Interaction  and 

Scattering  of  Euler  Ecuation  V-states  via  Contour  Dynamics,  Technical 
Report  ICMA-31-29,  Institute  for  Computational  Mathematics  and 
Applications,  Department  of  Mathematics  and  Statistics,  University  of 
Pittsburgh,  Pittsburgh,  PA,  15261. 


49 


Picone,  J.M.  and  J.P.  Boris  (1983),  Vorticity  Generation  by  Asymmetric  Energy 
Deposition  in  a  Gaseous  Medium,  Phys .  Fluids  26(2),  365-382. 

Roshko,  A.,  (1976),  Structure  of  Turbulent  Shear  Flows:  A  Mew  Look,  AIAA 

Journal  14,  1349-1357. 

Turkel,  E.  (i980),  Numerical  Methods  for  Large-Scale  Time -Dependent  Partial 
Differential  Equations,  in  Computational  Fluid  Dynamics,  pp. 128-262  ed. 
W.  Kollman,  Hemisphere  Press,  Wash.,  D.C. 

Wygnanski,  I.,  and  H.  Fiedler  (1970),  The  Two-Dimensional  Mixing  Reqion, 

J.  Fluid  Mech.  41,  327-361. 

Yule,  A.J  (1981)  Investigations  of  Eddy  Coherence  in  Jet  Flows,  in  The  Role 
of  Coherent  Structures  in  Modelling  Turbulence  and  Mixing,  ed . 

J.  Jiminez,  Sprinaer-Ver lag,  New  York. 

Zabusky,  N.J.  (1981),  Computational  Synergetics  and  Mathematical  Innovation, 
J.  Comp.  Phys.  43,  195. 


i 


